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WAVES IN ZND 



JEFFREY HUMPHERYS AND KEVIN ZUMBRUN 

Abstract. As described in the classic works of Lee-Stewart and Short-Stewart, the numerical 
evaluation of linear stability of planar detonation waves is a computationally intensive problem 
of considerable interest in applications. Reexamining this problem from a modern numerical 
Evans function point of view, we derive a new algorithm for their stability analysis, related to 
a much older method of Erpenbeck, that, while equally simple and easy to implement as the 
standard method introduced by Lee-Stewart, appears to be potentially faster and more stable. 



1. Introduction 

As described for example in [201 ED E2J [Ml EH H5] , the numerical stability analysis of detona- 
tion wave solutions of the Zeldovich-von Neumann-Doring (ZND), or reactive Euler equations, 
is a rich and computationally challenging problem. Planar detonation waves can often change 
stability as physical parameters are varied, undergoing interesting bifurcations to pulsating, 
spinning, and cellular solutions [IH [23l El [321 E3 [291 SZl [13 H9] . This motivates the numerical 
study of their stability, originated by Erpenbeck in |2Q(. I2XJ . both for its interest in its own right 
and as a benchmark for more general time-evolution codes \12\ H5] . 

Due both to the number of physical parameters (four for a polytropic gas3) and the difficulty of 
individual computations, this problem has proven to be numerically intensive. In their classical 
1990 paper [M], in which they introduced the algorithm that has become the modern-day 
standard, computing accurately for the first time the stability boundaries for one-dimensional 
detonations, Lee and Stewart conclude (p. 131 of the reference): "Finally, we point out that 
though our scheme is direct and easy to implement, complete investigation of the various regions 
of parameter space is computationally intensive. Any equivalent or more efficient numerical 
method should be considered a valuable contribution and such approaches are needed to further 
explore the parameter regimes of instability." 

Despite these comments, the basic algorithm introduced by Lee-Stewart (or perhaps vari- 
ants thereof) as described in the 2006 survey [45] appears still to be the current state of the 
art. Of course, computational power has increased tremendously in the interim, making once- 
prohibitive computations now accessible. Nonetheless, it seems of interest to explore more 
efficient algorithms if they can be found. 
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In particular, the computations of [M] were carried out in 1990 on a Cray X-MP/48 super- 
computerU with several hours required to produce individual figures. (For example, Fig. 9 of 
|34j tracking the top 6 unstable eigenvalues of detonations of a polytropic gas with gas constant 
7 = 1.2 as activation energy is varied was reported to require 5 hours of computation.) Today, 
substantially more computing power is available in a standard desktop PC, and a relatively 
inexpensive multi-core workstation offers substantially moreJl Hence, the challenge is trans- 
posed from the level of the national lab to the level of individual users, and from feasibility to 
practical ease of use. However, the impetus is no less real to reduce computation time from 
hours to the minutes required for interactive numerical explorations, and such improvement 
would undoubtedly lead to further advances in our understanding of detonation phenomena. 

Meanwhile, in parallel development, there has been considerable activity, centered around 
the Evans function [U [39[ [25] . in the numerical evaluation of stability of viscous shock waves 
and other traveling front or pulse and boundary layer solutions arising in a variety of equations 
[HI US [HI H3 EH HI EH Ell El El El [HI IS], some of which problems- see, e.g., [271 H] 
exhibit complexity rivalling that of detonations. The authors and collaborators have developed 
a general model-independent method and set of numerical principles for the treatment of such 
problems |28[ [54"] . encoded in the MATLAB-based platform STABLAB [6], which performs 
extremely well on all of the above-described applications. 

At the same time, there has been a successful push to place detonation stability in a common 
framework with stability of shock waves [50l [35l EH [291 SSI E21 ES] . In particular, in [50l 1291 
l4"9 j I52 [ [55]. the determination of stability of both viscous (reactive Navier-Stokes) and inviscid 
(reactive Euler or ZND) detonations has been reduced to the computation of an Evans function 
defined exactly as in the viscous shock and other cases described above. Thus, it is a natural 
step to study ZND stability within this common framework, using the general tools of |28[ 154], 

In this paper, we do exactly that, proposing a new algorithm for the numerical determination 
of stability of ZND detonations derived from the point of view of [28j [54] • Surprisingly, though 
both are shooting methods, this is quite different from the Lee-Stewart algorithm currently in 
standard use, shooting from x = — oo to x = rather than from x = to x = — oo as in 
[34] : indeed, it is more closely related to the original algorithm of Erpenbeck [21] . The precise 
relations between the various methods are described in Section [4] 

The advantage of shooting from — oo to is that we seek generalized eigenfunctions decaying 
exponentially at — oo. Thus, in the forward direction (— oo — > 0), the desired solution grows 
exponentially, while error modes are exponentially damped. By contrast, integrating in the 
backward direction (0 — > — oo), the desired solution decays exponentially while error modes 
are exponentially amplified, a numerically undesirable situation ("numerical pitfall 1" of [54J). 
For this reason, we expect that our algorithm should be faster and better conditioned than the 
Lee-Stewart algorithm currently in use. However, there are other aspects that cloud the issue, 
in particular the singular perturbation structure that arises in the high-activation energy or 
"square- wave" limit in which instabilities are often studied [20 [ I22 [ [23 j I17 [ [2]. For this reason, 



2 A Cray X-MP/48 cost roughly $15-20M dollars in the mid-1980's, having 2 processors with a 105 MHz clock 
speed and a theoretical peak performance about 200 MFLOPS per processor or 400 MFLOPS total. 

2010 Mac Pro 8-core (2 quad-core Xeon processors) for example is a $4-5k system with a 2.5GHz clock 
speed and a theoretical peak performance around 10 GFLOPS per core or 80 GFLOPS total. Hence, it has 
roughly 200 times the processing power at a five thousandth the price (not even adjusting for inflation). 
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careful comparison of methods in physically relevant regimes is an important step before making 
conclusions. 

In the present paper, we introduce the algorithm, and give some supporting numerical exper- 
iments for a simple model equation indicating the advantages of our approach. Followup work 
in [10] flT] indicates that, also in physically realistic settings, the algorithm performs favorably 
compared to the current standard. Specifically, the standard adaptive- mesh version of the al- 
gorithm described here appears to outperform the fixed-mesh algorithm described in |34] H5] by 
2-3 orders of magnitude. Much of this improvement appears to be due to the difference between 
fixed and adaptive mesh. However, even compared to an adaptive-mesh version of the method 
of Lee-Stewart, our algorithm appears to be 1-10 times faster, depending on the parameter 
regime: at the least, it is equivalent, and in some situations substantially more efficient. 

Plan of the paper. In Section[2j we review the ZND equations and detonation structure. In 
Section [3j we give a simple derivation of the Evans /Lopatinski function condition for detonation 
stability from a general point of view following [50, 29]. For clarity, we specialize in most of the 
discussion to the single-species, ideal gas case with Arrhenius ignition dynamics, working in the 
same framework as in |34j . The general case is discussed briefly in Remark 15.11 In Section [H 
we determine the relation between the derived Evans/Lopatinski condition the related stability 
determinants of Erpenbeck [20] and Lee-Stewart [34]. In Section [U we describe a proposed 
numerical implementation within the standard STABLAB package developed by the authors 
and collaborators. Finally, in Section [HJ we present numerical experiments for a simple model 
indicating the advantages of integrating in the forward direction and factoring out expected 
decay at — oo as prescribed in [2"5| [54"]. 

2. ZND DETONATIONS 

2.1. The model. In Eulerian coordinates the Zeldovich-von Neumann-Doring (ZND) equa- 
tions of reacting gas dynamics in one space dimension may be written as 

Pt + {pu) x = 

{pu) t + (pu 2 +p) x = 

( P E) t + ((pE+p)u) x = 

(pY) t + (puY) x = -p<p(T)KY, 

where p, u, p, E, T G R 1 represent density, velocity, pressure, total energy, and temperature, 
and Y = (Y%, . . . , Y r ) € W the mass fractions of reactantsQ Here, E = v? /2 + e is the non- 
reacting gas-dynamical energy E = u 2 /2 + e modified by chemical potential according to 

e = e + qY, 

where e is the specific internal energy of the gas and qY is the specific chemical energy. The 
matrix K € W xr and vector q € M} xr measure the rates of reaction and the heat released in 
reaction, respectively, and (p is an "ignition function" that is positive for T above some ignition 
temperature Tj and zero for T < Tj, serving to "turn on" the reaction. The matrix —K is 
assumed to be stable, i.e., to have spectrum of strictly negative real part, so that reaction in a 



Alternatively, the equations may be written in terms of progress variables Aj ; = 1 — Yj 22, 34l 136] . 
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quiescent flow indeed proceeds to the completely burned state Y = 0. In the simplest case of a 
single-species, exothermic reaction, Y € M 1 is a scalar, and K and q are positive constants. 

The system is closed by specifying equations of state (i.e., thermodynamic relations) p = 
p(p, e, Y) and T = T(p, e, 1") and the ignition function. Standard assumptions (in particular, 
the ones made in [34J, etc.) are the ideal gas laws 



where T, C v > are constants determined by the nature of the gas, and the modified Arrhenius 
law 



where Ea is the activation energy, R = 7CV is the gas constant, and (3 is an artificial smooth 
cutoff function with the property that (3 = 1 for T >T l and /3 = for T < TjJ^ Under usual 
assumptions, the specific form of the function /3 plays no role in the analysis; see Remark 12. 21 

Remark 2.1. More realistic rate laws r(p,T,Y) may be considered in place of the linear law 
r = —pip(T)KY with little additional difficulty |34j ; however, we lose the explicit form of the 
reaction profile (|5.3p computed in Section \5.1\ In the single- species case, these are equivalent. 

2.2. Alternative formulation. Subtracting q times the fourth equation of (12. ip from the 
third equation, we obtain the alternative formulation 



in terms of the usual gas-dynamical variables p, u, E. We alternate between the two formula- 
tions as convenient for the analysis. 

2.3. Detonation waves. For temperatures T < Ti below igition level, equations (12.ip evi- 
dently reduces to the usual Euler equations of nonreactive gas dynamics, with the reactants Y 
convected passively by the velocity field u. In particular, so long as T(p±, e±, Yq) < Ti, they 
support as traveling-wave solutions ordinary gas-dynamical shock waves 



(2.2) 




(2.3) 




(2.4) 



Pt + {pu) x = 
(pu) t + (pu 2 +p) x = 
(pE) t + ((pE+p)u) x = pq^(T)KY 
(pY) t + (puY) x = -p<p(T)KY 




x - st > 



satisfying the Rankine-Hugoniot conditions 






or, equivalent ly, 



(2.6) 




2 





The latter, standard modification circumvents the "cold-boundary difficulty" that the unburned state Y = 1 
is not an equilibrium for the exact Arrhenius law (3 = 1, and so steady traveling detonation waves do not exist. 
Though not mentioned, this assumption is also made implicitly in [34], etc. 
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where for an arbitrary function h(p,u, E,Y), [h] := h(p + ,u + , E + ,Y + ) — h(p_,u_, E_,Y_) de- 
notes jump across the discontinuity. This also holds if there is no reactant, Yo = (0, . . . , 0). 

If, on the other hand, Yq ^ (0, . . . , 0), and T + < T, but T_ > Tj, with u± < s (alternatively, 
T + > Ti and T_ < Tj, with u± > s), then there appears a different type of traveling- wave 
solution known as a strong detonation, given by (z = x — st) 



(2.7) (p,u,E,Y)(z) 



j (p + ,u + ,E + ,Y ) z>0 
\(p,u,E,Y)(z) z<0, 

where Y(z) satisfies the smooth traveling-profile ODE 

(2.8) m-s)Y)' = -(xp{T)KY 

on ( — oo , 0], with initial condition Y(0) = Yq, decaying to the completely burned state (0, . . . , 0) 
as z — > — oo, with (p,u,E) = (p,u, E)(Y) determined through the generalized Rankine- 
Hugoniot relations 

sp — pu = (sp — pu)± 

(2.9) spu - (pu 2 + p) = (spu - (pu 2 + p))± 

spE — (pE + p)u = (spE — (pE + p)u)± 

obtained by integrating the remaining traveling-profile equations 

(sp — pu)' = 

(2.10) (spu - (pu 2 +p))' = 

(spE - (pE + p)u)' = 

from to z (where z < 0) and recalling the Rankine-Hugoniot conditions (|2.6p satisfied across 
the jump at z = 0. 

That is, strong detonations moving to the right with respect to fluid velocity u (i.e., u < s, 
where s is the speed of the detonation) have the structure of an initiating gas-dynamical shock 
called the Neumann shock, which rapidly compresses the gas, raising temperature to the point 
of ignition, followed by a reaction zone (the profile (p,u, E,Y)) resolving to the final burned 
state. This characteristic "detonation spike" in temperature and pressure profiles agrees well 
with observed features in laboratory experiments. 

Substituting into fj2.8f) the first relation in (12. 9p and introducing the constant m := (p(s — 
u ))±i we obtain the simplified reaction equation 

(2.11) Y' = mr x pv(T)KY 

that we will actually use to solve for the profile. Further simplifying (I2.9p . we obtain 

sp — pu = (sp — pu)± 

(2.12) spu — (pu 2 + p) = (spu — (pu 2 + p))± 

spE — (pE + p)u + mqY = (spE — (pE + p)u + mqY)±. 



An application of the Implicit Function Theorem reveals that (|2.9p (as, likewise, the original 
ODE (|2.10p ) may be solved for (p, u, E) in terms of Y so long as the gas-dynamical state (p, u, E) 
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remains noncharacteristic with respect to speed s, or, equivalently, the Rankine-Hugoniot re- 
lation (|2.12p remains full rank in (p,u,E). For typical reactions and equations of state, in 
particular ideal gas dynamics with single exothermic reaction, this condition holds for all so- 
lutions of (|2.12p with Yj > 0, except for special limiting values of s for which the asymptotic 
state Y = is characteristic, or "sonic"; see, e.g., [35j . These limiting, characteristic waves 
are called Chapman- J ouget detonations, and have a special place in the theory. The usual, 
noncharacteristic type are called overdriven detonations. 

For our present purposes, the main import of characteristicity is that the eigenvalue equation 
becomes singular at x — > — oo in the coordinates we use here, complicating the discussion. For 
simplicity, we restrict hereafter to the overdriven case. The Chapman- Jouget case may be 
treated similarly using ideas of [M]; see Remark 15. II 

Remark 2.2. For the modified Arrhenius ignition function (|2.3I) . a standard assumption is 
that T > T l , all x < 0, T + < Tj, so that (3 = 1 for x < and /3 = for x > 0. Under this 
assumption, the specific form of the cutoff f3 plays no role in the analysis. 



3. Linearized stability analysis: the Evans-Lopatinski determinant 



We now carry out a linearized interface analysis, loosely following [29jjj Setting V := 
{p,u,e) T , write (|2.4j) in abstract form as 

(3.1) F°(W) t + F 1 (W) x = R(W), 
W, F*, Re R 3+r , where 

(3.2) W.= (V), F>:=(m), « :- 



Y)> r ■ \Yg?{V))> n - y-KYiP(W)J 

pu 

/ U :=l pu I, / x :=l pu 2 +p(p,e,Y) 
(3.3) \p(e + n 2 /2)/ \(p(e + u 2 /2) + p(p, e, Y))u, 

,o. „ i ^ ... 



{ " \ 




PU 


\p(e + U y2) 





9 ■= P, 9 = pu, Q := ( , ip := p4>(T(p,e,Y)). 



with V, ft G M 3 , Y € W, g* , ip E R 1 , Q E 



ixr 



Remark 3.1. A minor departure from \19\ [20| [50| [29] is to admit the possible dependence of 
pressure and temperature on chemical makeup of the gas (Y) , an important feature in realistic 
modeling of reactive flow. 

To investigate solutions in the vicinity of a discontinuous detonation profile, we postulate 
existence of a single shock discontinuity at location X(t), and reduce to a fixed-boundary 
problem by the change of variables x — >• x — X(t). In the new coordinates, the problem 
becomes 

(3.4) F°(W) t + (F 1 (W)-X'(t)F°(W)) x = R(W), x + 0, 



^See also the related [501 1351 [36] . and the original treatments in [191 1201 [34] . etc. 
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with jump condition 

(3.5) X'(t)[F°(W)} - [F\W)} = 0, 

[h(x,t)] := h(0 + ,t) — h(0~,t) as usual denoting jump across the discontinuity at x = 0. 

3.1. Linearization. Without loss of generality, suppose for simplicity that the background 
profile If is a steady detonation, i.e., s = 0, hence (W,X) = (W, 0) is also a steady solution of 
(|3.4p ~ (|3.5p . Linearizing (|3.4p - (|3.5p about the solution (W, 0), we obtain the linearized equations 

(3.6) A°{W t - X'{t)W'(x)) + (A 1 W) X = CW, 

(3.7) X'{t)[F°(W)]-[A 1 W} = 0, x = 0, 

(3.8) A j := (d/dW)Fi, C := (d/dW)R. 



3.2. Reduction to homogeneous form. As pointed out in [29j, it is convenient for the 
stability analysis to eliminate the front from the interior equation (|3.6p . Therefore, we reverse 
the original transformation to linear order by the change of dependent variables 

(3.9) W->W-X(t)W'(x), 
following the calculation 

W(x - X(t),t)) - W{x,t) ~ X(t)W x {x,t) ~ X(t)W'{x). 

approximating to linear order the original, nonlinear transformation. Substituting (|3.9p in 
(I3.6p - (l3.7p . and noting that x-differentiation of the steady profile equation F 1 (W) X = R(W) 
gives 

(3.10) (A 1 (W)W'(x)) x = C(W)W'(x), 
we obtain modified, homogeneous interior equations 

(3.11) A°W t + (A 1 ^. = CW 

agreeing with those that would be obtained by a naive calculation without consideration of the 
front, together with the modified jump condition 

(3.12) X'(t)[F°{W)} - X{t)[A l W'{x)\ - [A^] = 

correctly accounting for front dynamics. 

The reduction to homogeneous interior equations puts the linearized problem in a standard 
linear boundary- value-problem format for which stability may be investigated in straightforward 
fashion by the construction of an Evans/Lopatinski determinant. Besides simplifying consider- 
ably Erpenbeck's original derivation of his equivalent stability function [20], the homogeneous 
format makes possible the application of standard numerical Evans function techniques for its 
evaluation. This useful reduction was first carried out, in slightly different form, in |29j . The 
transformation (I3.9P is of general use in interface problems, comprising the "good unknown" of 
Alinhac [3]. A similar discussion in the simpler context of shock waves may be found in |24j : 
however, in this case, W'(x) = 0, and so the transformation (|3.9p does not make itself evident, 
nor do front dynamics modify (j3. 12|) . 
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3.3. The stability determinant. Seeking normal mode solutions W(x,t) = e M W(x), X(t) = 
e xt X, W bounded, of the linearized equations (|3.1ip ^ (|3.12p . we are led to the generalized 
eigenvalue equations 

(A 1 W)' = (-\A° + C)W, x^O, 

X(X[F°{W)} - [A^'ix)]) - [A X W] = 0, 
where "/" denotes d/dx, or, setting Z := A X W, to 

(3.13) Z' = GZ, x + 0, 

(3.14) X(X[F°(W)} - [A l W'{x)}) - [Z] = 0, 
with 

(3.15) G:= (-AA^C)^ 1 )- 1 . 

Here, we are implicitly using the following elementary observation. 

Lemma 3.2. A l (W(x)) is invertible for all x such that df/dV is invertible (i.e. V is non- 
characteristic as a gas- dynamical state with Y held fixed). 

Proof. Similarly as in the discussion of existence of steady profiles, we may by subtracting Y 
times the first row of A 1 from the block Y-row, reduce A 1 to block upper-triangular form, with 
diagonal blocks df/dV and g^V, Y)I rxr with g^V, Y) = pu / 0. □ 

Remark 3.3. As discussed in Section \2.1\ this assumption is essentially necessary already for 
existence of a steady profile. In particular, it is satisfied for the usual ideal gas equation of state. 

We require also the following fundamental properties. 

Lemma 3.4 ( |19^ [20| I29j). On 5ReA > 0, the limiting (3 + r) x (3 + r) coefficient matrices 
G± := lim^ioo G(z) have unstable subspaces of fixed rank: full rank 3 + r for G + and rank 
2 + r for G-. Moreover, these subspaces have continuous limits as KeA — > 0. 

Proof. Straightforward calculation using the fact that G± are block upper-triangular in (V, Y); 
see, e.g., [19], [20j [50l El] in the case that /, g depend only on V . □ 

Corollary 3.5 ([50j[29]). On KeA > 0, the only bounded solution of (|3. 13[) for x > is the triv- 
ial solution W = 0. For x < 0, the bounded solutions consist of an (r + 2) -dimensional subspace 
Span {Z^ , . . . , Z^ +2 }{X, x) of exponentially decaying solutions, analytic in A and tangent as 
x — > —oo to the subspace of exponentially decaying solutions of the limiting, constant- coefficient 
equations Z' = G-Z; moreover, this subspace has a continuous limit as KeA — > 0. 

Proof. The first observation is immediate, using the fact that G is constant for x > 0. The 
second follows from asymptotic ODE theory, using the "gap" or "conjugation" lemmas of 
[25"1 [30] , [38] together with the fact that G decays exponentially to its end state as x — > — oo. 
See [291 [521 [55] for details. □ 
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Definition 3.6. We define the Evans-Lopatinski determinant 

L>(A):=det(^-(A,0), Z- +2 (A,0), \[F°(W)} - [A l W'(x)}) 

= det(Z 1 -(A,0), Z f 7 +2 (A,0), X[F°(W)] + ^'((T)) , 

where Z^(X,x) are as in Corollary 13,51 

The function -D is exactly the stability function derived in a different form by Erpenbeck |20j : 
see Section 14.21 below. The formulation (|3.16p is of the standard form arising in the simpler 
context of (nonreactive) shock stability [37} I19j. Evidently (by (|3.14p combined with Corollary 
I3.5p . A is a generalized eigenvalue/normal mode for SfteA > if and only if D(X) = 0. 



Remark 3.7. As noted in [521 155] . consideration of the traveling-wave equation F(W)' = 
AW' = R{W) yields the simpler formula 

(3.17) D(X) = det(Z 1 -(A,0), Z p - +2 (A,0), \{F°(W)] + R(W(0~))) . 

3.4. Dual formulation. The (n + r) x (n + r) determinant (I3.17P may be expressed more 
succinctly in dual form 

(3.18) D(X) = Z~(X, 0) • {\[F°(W)\ + R(W(0~)), 
where Z~(X, x) is the cross product Z^ A • • • A Z r+2 (A, x) defined by 

Z~ ■ x = det (Zf, ■ ■ Z~ +2 , x) . 

The vector Z~ may alternatively be characterized directly as the unique up to constant factor 
bounded solution on x < of the adjoint ODE 

(3.19) Z' = -G*Z, 

which, as x — > — oo is both exponentially decaying and tangent to the corresponding expo- 
nentially decaying one-dimensional subspace of bounded solutions of the limiting constant- 
coefficient equations Z' = —G*_Z. It may be specified analytically in A by the additional 
requirement 

(3.20) Il(Z-) conj (-M) = £(X), 

M > 0, where I is an analytically chosen left eigenvector of G_(A) associated with the unique 
eigenvalue g~(X) of negative real part and IT the associated eigenprojection. Here, and else- 
where, con3 denotes complex conjugate. By (I3.20P together with the tangency property, Z~ is 
well-approximated at x = —M, for M > sufficiently large, by 

(3.21) Z-(-M) =£ conj {X). 

This reduces the approximate evaluation of D(-) to the straightforward and extremely well- 
conditioned numerical problem of integrating a single exponentially growing (in forward direc- 
tion) mode from x = —M to x = 0. The stability of the computation derives from the fact that 
errors lying in other, exponentially decaying modes, are exponentially damped |54j . 

Alternate initialization. Alternatively, following \14 \ \15 \ [T6 | 113]. Z~ may be specified by 
boundary conditions at — oo, via 

(3.22) lim e 9 -"' x Z-{x) =e conj {X), 
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whence (|3.2ip becomes 

(3.23) Z-(-M) = e - g -° n3M e conj (X). 

This is the method that we prescribe here. It has the advantage of removing the dependence 
of Z~ on the artificial parameter M, allowing the flexible choice of M in different parameter 
regimes, as dictated by numerical considerations, while preserving analyticity. However, in 
practice, there is usually not much difference between (|3,21|) and (|3.23p . In particular, if, as in 
[34] . one is not interested in analyticity, then one may vary M freely in (13.211) as well. 

4. Relations to other methods 

4.1. Relation to the method of Lee and Stewart. Denoting by Zq the solution on x < 
of the forward eigenvalue ODE (pTT3j) with initial conditions Z (0) := X[F°(W)} + R(W(0~)), 
we have by standard duality properties that 

(4.1) Z- ■Z (X,x) = D(X) 

is independent of x < 0, or {Z~ ■ Zq)'(X, x) = 0. Taking x = —M and recalling (|3.2ip . we arrive 
at the alternative Evans-Lopatinski approximation 

(4.2) D(X) ~ £ con i(X) ■ Z Q (X, -M) 

used by Lee and Stewart [34] , where i con ^ ■ Zq(—M) = is their "nonradiative condition" 
enforcing boundedness of Zq. The solution of Zq from x = to x = —M, on the other hand, is 
numerically comparatively ill-conditioned in the vicinity of roots of -D(-), since Zq in this regime 
is approximately exponentially decaying in the backward direction while errors are exponentially 
growing^] The version (13.18P is therefore much preferable from the numerical point of view, at 
least when used (as here, and in [34]) as a shooting method. 

4.2. Relation to the method of Erpenbeck. Erpenbeck [21] computes Z~ in much the 
same way as we do here. However, in place of the homogeneous duality relation (14. ip . he uses 
the "inhomogeous Abel relation" 

(4.3) (Z- ■Zq)'(X,x)) = Z XW'(x), 

valid for the solution Zq of the inhomogeneous equation Z' = GZ + XW'{x) with initial data 
Zo(0) := X[W] deriving from the unmodified equations f)3.6f) ()3. 7j) . together with W'(-oo) = 0, 
to evaluate 

D(X) = / Z~{y) ■ XW\y)dy + Z"(0) • X[W\. 

J — oo 

Though it is mathematically equivalent to the homogeneous scheme described above, this has 
the disadvantage that it is difficult to implement adaptive control on truncation error simul- 
taneously for the ODE and quadrature steps. Indeed, the method is in general a bit more 
cumbersome to implement and understand than either of the previous two described methods. 
As a one-time cost, the latter is a rather minor point. However, the implications of the for- 
mer for performance appear to be significant. Our experience in similar Evans function-type 

More precisely, they solve the inhomogeneous equations Zq — GZo + XW' (x) with initial data Zo(0) := A[W], 
and compute £ conj (X) ■ Zo(\, —M) ~ i coni> (A) • Zo(A, — M), which is numerically equivalent. Here we are using 
Zq — Zq — W'(x) — > as x — >• — oo. 
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shooting computations [T6j [2HJ [31 [26] of spectra of asymptotically constant-coefficient opera- 
tors is that a fixed-step scheme can be orders of magnitude slower than a comparable adaptive 
scheme; see [53] for a general discussion of performance of numerical Evans/Lopatinski solvers. 
Moreover, even in the solution of Z alone, the use of an adaptive solver without factoring out 
expected decay is much less effective in our experience ("numerical pitfall 3" of [54] ). 

4.3. Expression as boundary-value solver. We mention in passing an alternative "local 
Evans function" formulation in the spirit of [34j . suggested by Sandstede [41] as a general 
method for numerical Evans function investigations using collocation / 'continuation rather than 
shooting. By the analysis of the previous subsections, we may recast the eigenvalue equation 
(I3.13p - (13.14p as in [33] as an overdetermined two-point boundary-value problem Z' = GZ with 
r + 4 boundary conditions 



Relaxing at random one of the r + 3 conditions at x = 0, say the requirement on the jth 
coordinate, we generically obtain a well-posed boundary- value problem with the correct number 
r + 3 of boundary conditions; one of the coordinates will always suffice. More, the projective 
boundary-condition at x = — oo is numerically "correct", making this problem extremely well- 
conditioned for solution by collocation/continuation methods (see, e.g., [30])- Defining Z(X,x) 
to be the solution of this relaxed problem, we may then define a local, analytic Evans function 



that is numerically well-conditioned and vanishes if and only if A is an eigenvalue. This gives a 
second way to convert (I4.4h into a numerically well-conditioned problem, though the speed and 
simplicity of shooting is lost in this approach, along with global analyticity useful for winding 
number calculations. We shall not investigate this method here, but note that it could be useful 
in extreme conditions such as the ultra-high activation energy limit [17] . 



We now describe in detail the numerical algorithm proposed to compute (|3. 18[) . following the 
general approach set out in [TBJ [28j [531 [53] • 

5.1. Computing the profile. In Evans function computations, a delicate aspect is often the 
computation of the background nonlinear profile. We sidestep this issue by the explicit so- 
lution technique used in [191 [2Q1 E3]) modified slightly to accomodate the multi-species case 
(specifically, the simplified uniform ignition one considered here). 
Introducing the new variable y defined by 



(4.4) 



Z(0) := X[F°(W)\ + R(W(Q-)) 



lim l conj ■ Z{x) = 0. 



D(X) := ej ■ (Z(X, 0) - (X[F°(W)} + R(W(0~)))) 



5. Numerical implementation 




±, we reduce the reaction equation (|2.1ip to 
dY/dy = KY, Y(0) = Y , 



dy/dx = m p<p(T), y(0) = 0, 



Y(y) = e K yy 
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from which the full profile can be recovered through (|2.12p . either by explicit calculation, as 
carried out for ideal gas dynamics in Appendix [Bj or, more generally, by Newton iteration. 

Remark 5.1. In the single- species case, (|5.2p reduces to the change of coordinates x — > y := 
logy used in [34] ; general, nonlinear rate laws, or Chapman- J ouget waves, may be accomodated 
by a change of variables x — > Y r for appropriate r, as discussed in |34j . 

5.2. Computing the stability determinant. The linearized stability analysis can then be 
carried out in the variable y defined in (|5.ip , using the instantaneous change of variables formula 

(5.4) dx/dy = m/(pip(f )), y < 0. 

Remark 5.2. Since the righthand side of (15, 4p is uniformly positive and bounded, the variables 
x and y are equivalent in the sense that Cx < y < x/C for x, y < 0, for some C > 0. 

Specifically, we solve from y = -M to y = the ODE (d/dy)Z = -m/(p<p(T))G*(y)Z, with 

initial condition Z~(-M) = e' 9 ^ WM £ conj (X), M > sufficiently large, where the vector £(X) 
and limiting eigenvalue g\ = +c_) -1 are as computed in eqs. (|A.2p and (|A.3P of Appendix 
|A"1 the coefficient G(A, V, Y) is as described in eqs. (|3. 15|) . (j3.8|) . and ()3.2j) ()3.3|) . and the profile 
(V,Y)(y) is as computed in Appendix [Bj As prescribed in (|3.18p . we may then compute the 
stability determinant D(X) = Z~(\,0) ■ (X[F°(W)] + R(W'(Q-))). 

More precisely, we may solve the numerically more advantageous equations 

(5.5) (d/dy)Z = -(m/ptp(f))(G(y) + 9 -(X)I)*Z, 

with initial conditions Z{y) := e~^ mg ~ ^ v ^ y Z(y), and compute 

D(X) = Z~(X, 0) • (X[F°(W)] + R(W(0-))). 

This may readily be computed with good results by an adaptive solver such as the standard 
RK45; see [HI [281 El] for further discussion. 

5.3. Determination of stability: winding number vs. stability curves. With an Evans 
solver in hand, stability may be checked either by winding number computations as in |214l4"l[26]. 
or by root-following methods based on the Implicit Function Theorem, as in [34J. In the first 
method, a large semicircle S centered at the origin and lying in 5RA > is mapped by D, and 
the number of zeros of D (unstable normal modes) lying within S computed using the principle 
of the argument, making use of the underlying analyticity of D. Unstable modes lying outside 
S may be excluded by a separate, asymptotic, argument based on high-frequency behavior of 
D p31 [151 [2S]; for implementations in the context of ZND, see [551 E3] (analytical) or [101 [XT] 
(numerical). In the second method, individual roots are followed, avoiding the need to compute 
around a contour, but typically requiring an extra Newton iteration with each change in model 
parameters; see, for example, |34t 144] . Both are by now completely standard. 

6. A SIMPLE MODEL PROBLEM 

We conclude by an examination of efficiency within the context of a simple but illustrative 
model problem. Consider the ODE 

(6.1) y' = A(x,X)y, A(x,X) = X^ 2x ^) 
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defined on — oo < x < 0, a; € R, A € C, yd C 2 , with boundary conditions y ~ e AX//2 (l, 0) T 
as x — > — oo and y(0) = (1,0) T , modeling a variable-coefficient eigenvalue problem of the form 
arising in ZND, where the coefficient c ^ encodes rapidity of exponential decay. As for 
ZND, the coefficient matrix is exponentially asymptotically constant as x — > — oo, with size 
growing linearly in A, and has a unique decaying mode as x — > — oo for all KA > 0, extending 
continuously to 5RA = 0. Thus, we may expect somewhat similar behavior, at least away from 
the high-activation energy "square-wave" regime. 

In this context, our proposed algorithm consists of factoring out the expected decay e Xx ^ 2 
from the solution to obtain a "neutral" equation 

(6.2) y' = A(x,X)y, A(x,X) = x(A B , 

y := ye Xx / 2 , then solving ()6.2[) from x = —M to x = and checking whether y(0) lies parallel 
to (1,0) T . For reasonable values of c, a computational domain of M = 5 is sufficient. The 
method of Lee-Stewart, consists roughly of integrating the original equation (16. ip from x = 
to x = — M; the method of Erpenbeck consists roughly of integrating (16. ip from x = —M to 
x = without first factoring out expected exponential decay. For comparison, we considered 
also a worst-case scenario with maximum amplification of error modes, integrating (|6.2p from 
x = to x = — oo. 

We computed all with the adaptive-mesh RK45 algorithm (ode45) supported in MATLAbJI 
with error tolerance set at the standard level 10 -5 used for Evans computations [26l EJ 
E], measuring efficiency by the number of mesh points /function calls required to complete 
the computation. Extreme cases are A real- the "best" case, with a spectral gap between 
exponentially growing and exponentially decaying modes at — oo- and A imaginary- the "worst" 
case from our standpoint, with neither spectral gap nor exponential decay. From the standpoint 
of the Lee-Stewart method, the best and worst cases would appear to be reversed. 

The results, displayed in Tables 1 and 2 for a typical value c = 10, indicate that the proposed 
new algorithm performs 1-5 times faster than (adaptive versions of) either the Erpenbeck or 
Lee-Stewart methods, depending on the value of A, with particular improvement as |A| becomes 
large. It should be noted, moreover, that this is only a comparison of speed (number of mesh 
points) for the various methods to produce output with fixed truncation error. If we consider 
also accuracy, i.e., convergence error, then the results could be expected to be more dramatic, 
since both Lee-Stewart and Erpenbeck methods are numerically less well-posed than the forward 
"neutral" algorithm that we propose. 



Appendix A. Calculation of £ 

In this appendix, we show how to calculate for general equations of state the initializing 
vector £(X) used in (I3.2ip . the unique stable left eigenvector of the limiting coefficient matrix 
(A.l) 

G- = (—XA V _ + C-)(A 1 _Y 1 



(-A<£-ifV-)(<7-) 



i \-i 



In practice, faster than corresponding fixed-mesh methods [541 1101 [TT] . 
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Table 1. Runs for Eq. (|6.2p . Forward corresponds to our proposed method, 
with expected decay factored out. Backward is a worst-case scenario not corre- 
sponding to any of the methods considered. 
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Table 2. Runs for Eq. (|6.ip . Forward corresponds to Erpenbeck method, 
backward to Lee-Stewart method. 



where for a general function h(V,Y), we use h- to denote h(V-,Y-). Here, we have strongly 
used Y— = to obtain the simple upper block-triangular form. 

By the upper block-triangular form of and the fact that the lower right-hand block 
has spectrum of positive real part for 9ieA > (since g° > always, g 1 < for right-moving 
detonations, and —K is assumed to have spectrum of negative real part), we find that £ T must 
be of form (iy,^), where ty is the unique unstable eigenvector, associated with eigenvalue a, 
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of the purely gas-dynamical matrix f v _(f v _) , and 

(A.2) £ Y = %(\$_(fv-)- 1 fY- + QK1>-)(\(g°_ - agl)I + K^)' 1 . 

To determine a, £y, and thereby £y, we observe that fy_{f v _)~ l , is related by similarity 
tranform M — > {fy_)~ l M f v _ to the inverse of the hyperbolic convection matrix 

{fy)~ l fy of the nonreactive Euler equations Vt + (fy) fyV x = written in nonconservative 
form in V coordinates with Y = 0. Thus, a -1 is an eigenvalue of (fy)^ 1 fy, be., a hyperbolic 
characteristic speed of the non-reactive Euler equations, and £y = where £ is the 

associated left characteristic direction (eigenmode). 

Noting that a -1 , as the unique positive characteristic at state V = V-, must be the largest 
characteristic speed, we have by standard formulae [HI 321 321 [351 EI] or direct calculation 

(A.3) £l = F(f v _)- 1 = (p p -cu + p- 1 p e (u 2 /2-e),c-p- 1 p e u,p- 1 Pe ), 

determining £ T (\) = (£y,£ Y ){^) through (|A.2|) . Note that £y is independent of A. For Y- 
independent equations of state, (|A.2p simplifies considerably, to £ Y = £yQKip_(\(g°^ —ag l _)I + 

Remark A.l. Noting that the e-component p~ 1 p e of £y does not vanish in (|A.3j) . we may 
alternatively rescale by p/p e to obtain an analytic choice of form £ T = (*, *, 1, *) convenient for 
numerical solution. 

A.l. Alternative, numerical computation. Alternatively, an analytic choice of £ may be 
determined numerically by solution of Kato's ODE [31] as described in [16 \ \28 \ [53l 154) . For 5?A 
bounded from zero, this involves finding numerically at each A-value the unique stable left and 
right eigenvectors of G_ and computing the associated eigenprojection for use in the Kato ODE 
as in the general problem-independent method of [16, 28, 53, 54|. At or near 3?A = 0, however, 
this method must be modified, since the stable eigenvector becomes neutral at 5RA = 0. A 
simple resolution is to notice that, there, the eigenvalues of G_ consist of a single eigenvalue 
with strictly positive real part, which may be discarded, and three eigenvalues of form gj = ctj\, 
where ctj (see above) are hyperbolic characteristic speeds for the non-reactive Euler equations, 
of which the one for which gj/X = aj < is the one associated with £. 



Appendix B. Ideal gas profile 

In this appendix, we explicitly solve (12. 12ft for the case of an ideal gas. Restricting to a 
steady shock, s = 0, and using the ideal gas law (12. 2j) . we may rewrite (12.12P as 

pu = p±u± := — m 

u + T- = u± + T — :=b 
(B.l) u u± 

U 2 - 1la_ 

Y + (r + i)e + 9 y = -^ + (r + i)e± + gy± :=c. 
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Combining the second two equations and simplifying gives (T+2)u 2 —2(T+l)bu+2T(c—qY) = 0. 
Solving using the quadratic formula, we obtain 



(B.2) 



m bu — u 2 



u 




where we have chosen the negative solution branch for u in accordance with the fact that [u] > 0, 
or, equivalently, [p] < 0, for a right-moving gas-dynamical shock, so that u(0~) < u + . (Recall 
that u(0~) and ii+ are the two branches of the square root for Y = 1, corresponding to the 
solutions of the Rankine-Hugoniot conditions for a nonreacting gas-dynamical shock.) With 
(|5.3p . (|B.2p gives an explicit expression for the profile as a function of variable y. 

For a given Neumann shock, there is a one-parameter family of possible endstates (p,u,e)- 
determined by the value of q, the maximum value of q corresponding to a Chapman- Jouget 
wave, for which the argument of the square root vanishes for y = 0. 
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